#include "../tree_fitch.h"
#include "phylip.h"
#include "dist.h"
#include <algorithm>
#include <stdexcept>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <math.h>
#include <matrix.h>
#include <ultoa.h>
using std::vector;
using std::string;
using std::istream;
using std::ostringstream;
using std::istringstream;

namespace fitch_tree
{
using namespace bcpp_matrix;

/* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#define zsmoothings     10    /* number of zero-branch correction iterations */
#define epsilonf        0.000001   /* a very small but not too small number  */
#define delta           0.0001      /* a not quite so small number */




t_fitch_params::t_fitch_params(	unsigned njumble_,
								eJUMBLE jumble_,
								eGLOBAL global_,
								eMIN_EV minev_,
								eNEGLEN negallowed_,
								double power_,
								unsigned outgrno_)
				:
			outgrno(outgrno_),
			njumble(njumble_),
			global(global_ == eDO_GLOBAL),
			jumble(jumble_ == eDO_JUMBLE),
			minev(minev_ == eDO_MIN_EV),
			negallowed(negallowed_ == eDO_NEGLEN),
			power(power_)
{
	if (minev)
		negallowed = false;
}






void secondtraverse(node *q, double y, long *nx, double *sum)
{
	/* from each of those places go back to all others */
	/* nx comes from firsttraverse */
	/* sum comes from evaluate via firsttraverse */
	double z=0.0, TEMP=0.0;

	z = y + q->v;
	if ( q->tip )
	{
		TEMP = q->d[(*nx) - 1] - z;
		*sum += q->w[(*nx) - 1] * (TEMP * TEMP);
	}
	else
	{
		secondtraverse(q->next->back, z, nx, sum);
		secondtraverse(q->next->next->back, z, nx,sum);
	}
}  /* secondtraverse */


void firsttraverse(tree& curtree, node *p, long *nx, double *sum)
{
	/* go through tree calculating branch lengths */
	if ( curtree.minev && (p != curtree.start) )
		*sum += p->v;
	if ( p->tip )
	{
		if ( !curtree.minev )
		{
			*nx = p->index;
			secondtraverse(p->back, 0.0, nx, sum);
		}
	}
	else
	{
		firsttraverse(curtree, p->next->back, nx,sum);
		firsttraverse(curtree, p->next->next->back, nx,sum);
	}
}  /* firsttraverse */


double evaluate(tree& curtree)
{
	++curtree.dots;
	double sum=0.0;
	long nx=0;
	/* evaluate likelihood of a tree */
	firsttraverse(curtree, curtree.start->back ,&nx, &sum);
	firsttraverse(curtree, curtree.start, &nx, &sum);
	curtree.likelihood = -sum;
	return(-sum);
}  /* evaluate */


void nudists(node *x, node *y)
{
	/* compute distance between an interior node and tips */
	long nq=0, nr=0, nx=0, ny=0;
	double dil=0, djl=0, wil=0, wjl=0, vi=0, vj=0;
	node *qprime, *rprime;

	qprime = x->next;
	rprime = qprime->next->back;
	qprime = qprime->back;
	ny = y->index;
	dil = qprime->d[ny - 1];
	djl = rprime->d[ny - 1];
	wil = qprime->w[ny - 1];
	wjl = rprime->w[ny - 1];
	vi = qprime->v;
	vj = rprime->v;
	x->w[ny - 1] = wil + wjl;

	// CHANGE HERE
	if ( wil + wjl < 0.0 )
		x->d[ny - 1] = 0.0;					// LEO ignore adjacent null connections
		//x->d[ny - 1] = 1000.0;
	else
		x->d[ny - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
	nx = x->index;
	nq = qprime->index;
	nr = rprime->index;
	dil = y->d[nq - 1];
	djl = y->d[nr - 1];
	wil = y->w[nq - 1];
	wjl = y->w[nr - 1];
	y->w[nx - 1] = wil + wjl;
	if ( wil + wjl < 0.0 )
		y->d[nx - 1] = 0.0;					// LEO ignore adjacent null connections
		//y->d[nx - 1] = 1000.0;
	else
		y->d[nx - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
}  /* nudists */


void makedists(node *p)
{
	/* compute distances among three neighbors of a node */
	long i=0, nr=0, ns=0;
	node *q, *r, *s;

	r = p->back;
	nr = r->index;
	for ( i = 1; i <= 3; i++ )
	{
		q = p->next;
		s = q->back;
		ns = s->index;
		//CHANGE HERE?
		if ( s->w[nr - 1] + r->w[ns - 1] < 0.0 )
			p->dist = 0.0;					// LEO ignore adjacent null connections
            //p->dist = 1000.0;
		else
			p->dist	= (s->w[nr - 1] * s->d[nr - 1] + r->w[ns - 1] * r->d[ns - 1]) /
					  (s->w[nr - 1] + r->w[ns - 1]);
		p = q;
		r = s;
		nr = ns;
	}
}  /* makedists */


void makebigv(node *p)
{
	/* make new branch length */
	long i=0;
	node *temp, *q, *r;

	q = p->next;
	r = q->next;
	for ( i = 1; i <= 3; i++ )
	{
		if ( p->iter )
		{
			p->v = (p->dist + r->dist - q->dist) / 2.0;
			p->back->v = p->v;
		}
		temp = p;
		p = q;
		q = r;
		r = temp;
	}
}  /* makebigv */


void correctv(tree& curtree, node *p)
{
	/* iterate branch lengths if some are to be zero */
	node *q, *r, *temp;
	long i=0, j=0, n=0, nq=0, nr=0, ntemp=0;
	double wq=0.0, wr=0.0;

	q = p->next;
	r = q->next;
	n = p->back->index;
	nq = q->back->index;
	nr = r->back->index;
	for ( i = 1; i <= zsmoothings; i++ )
	{
		for ( j = 1; j <= 3; j++ )
		{
			if ( p->iter )
			{
				wr = r->back->w[n - 1] + p->back->w[nr - 1];
				wq = q->back->w[n - 1] + p->back->w[nq - 1];
				if ( wr + wq <= 0.0 && !curtree.negallowed )
					p->v = 0.0;
				else
					p->v = ((p->dist - q->v) * wq + (r->dist - r->v) * wr) / (wr + wq);
				if ( p->v < 0 && !curtree.negallowed )
					p->v = 0.0;
				p->back->v = p->v;
			}
			temp = p;
			p = q;
			q = r;
			r = temp;
			ntemp = n;
			n = nq;
			nq = nr;
			nr = ntemp;
		}
	}
}  /* correctv */


void alter(node *x, node *y)
{
	/* traverse updating these views */
	nudists(x, y);
	if ( !y->tip )
	{
		alter(x, y->next->back);
		alter(x, y->next->next->back);
	}
}  /* alter */


void nuview(node *p)
{
	/* renew information about subtrees */
	long i=0;
	node *q, *r, *pprime, *temp;

	q = p->next;
	r = q->next;
	for ( i = 1; i <= 3; i++ )
	{
		temp = p;
		pprime = p->back;
		alter(p, pprime);
		p = q;
		q = r;
		r = temp;
	}
}  /* nuview */


void update(tree& curtree, node *p)
{
	/* update branch lengths around a node */

	if ( p->tip )
		return;
	makedists(p);
	if ( p->iter || p->next->iter || p->next->next->iter )
	{
		makebigv(p);
		correctv(curtree, p);
	}
	nuview(p);
}  /* update */


void smooth(tree& curtree, node *p)
{
	/* go through tree getting new branch lengths and views */
	if ( p->tip )
		return;
	update(curtree, p);
	smooth(curtree, p->next->back);
	smooth(curtree, p->next->next->back);
}  /* smooth */


void filltraverse(node *pb, node *qb, bool contin)
{
	if ( qb->tip )
		return;
	if ( contin )
	{
		filltraverse(pb, qb->next->back,contin);
		filltraverse(pb, qb->next->next->back,contin);
		nudists(qb, pb);
		return;
	}
	if ( !qb->next->back->tip )
		nudists(qb->next->back, pb);
	if ( !qb->next->next->back->tip )
		nudists(qb->next->next->back, pb);
}  /* filltraverse */


void fillin(node *pa, node *qa, bool contin)
{
	if ( !pa->tip )
	{
		fillin(pa->next->back, qa, contin);
		fillin(pa->next->next->back, qa, contin);
	}
	filltraverse(pa, qa, contin);
}  /* fillin */


void insert_(tree& curtree, node *p, node *q, bool contin_)
{
	/* put p and q together and iterate info. on resulting tree */
	double x=0.0, oldlike;
	hookup(p->next->next, q->back);
	hookup(p->next, q);
	x = q->v / 2.0;
	p->v = 0.0;
	p->back->v = 0.0;
	p->next->v = x;
	p->next->back->v = x;
	p->next->next->back->v = x;
	p->next->next->v = x;
	fillin(p->back, p, contin_);
	evaluate(curtree);
	do
	{
		oldlike = curtree.likelihood;
		smooth(curtree, p);
		smooth(curtree, p->back);
		evaluate(curtree);
	} while ( fabs(curtree.likelihood - oldlike) > delta );
}  /* insert_ */


void copynode(node *c, node *d)
{
	/* make a copy of a node */

	d->d = c->d;
	d->w = c->w;
	d->v = c->v;
	d->iter = c->iter;
	d->dist = c->dist;
	d->xcoord = c->xcoord;
	d->ycoord = c->ycoord;
	d->ymin = c->ymin;
	d->ymax = c->ymax;
}  /* copynode */


void copy_(tree *a, tree *b)
{
	/* make copy of a tree a to tree b */
	node *p, *q;
	b->spp = a->spp;
	b->nonodes = a->nonodes;
	for (unsigned i = 0; i < a->spp; i++ )
	{
		copynode(a->nodep[i], b->nodep[i]);
		if ( a->nodep[i]->back )
		{
			if ( a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1] )
				b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
			else if ( a->nodep[i]->back
					  == a->nodep[a->nodep[i]->back->index - 1]->next )
				b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
			else
				b->nodep[i]->back
				= b->nodep[a->nodep[i]->back->index - 1]->next->next;
		}
		else b->nodep[i]->back = NULL;
	}
	for (unsigned i = a->spp; i < a->nonodes; i++ )
	{
		p = a->nodep[i];
		q = b->nodep[i];
		for (unsigned j = 1; j <= 3; j++ )
		{
			copynode(p, q);
			if ( p->back )
			{
				if ( p->back == a->nodep[p->back->index - 1] )
					q->back = b->nodep[p->back->index - 1];
				else if ( p->back == a->nodep[p->back->index - 1]->next )
					q->back = b->nodep[p->back->index - 1]->next;
				else
					q->back	= b->nodep[p->back->index - 1]->next->next;
			}
			else
				q->back	= NULL;
			p = p->next;
			q = q->next;
		}
	}
	b->likelihood = a->likelihood;
	b->start = a->start;
}  /* copy_ */


void setuptipf(tree& curtree, const t_matrix_f& io, unsigned m, tree *t)
{
	/* initialize branch lengths and views in a tip */
	node *WITH;

	WITH = t->nodep[m - 1];
	std::copy(io.matrix[m - 1], io.matrix[m - 1] + t->nonodes, WITH->d.begin());
	for ( unsigned i = 0; i < t->spp; i++ )
	{
		if ( i + 1 != m && !io.is_missing(m-1, i))
		{
			if ( WITH->d[i] < epsilonf )
				WITH->d[i] = epsilonf;
			WITH->w[i] = 1.0 / exp(curtree.power * log(WITH->d[i]));
		}
		else
		{
			WITH->w[i] = 0.0;
			WITH->d[i] = 0.0;
		}
	}
	for (unsigned i = t->spp; i < t->nonodes; i++ )
	{
		WITH->w[i] = 1.0;
		WITH->d[i] = 0.0;
	}
	WITH->index = m;
	if ( WITH->iter ) WITH->v = 0.0;
}  /* setuptipf */


void buildnewtip(tree& curtree, const t_matrix_f& io, unsigned m, unsigned nextsp)
{
	/* initialize and hook up a new tip */
	node *p;
	setuptipf(curtree, io, m, &curtree);
	p = curtree.nodep[nextsp + curtree.spp - 3];
	hookup(curtree.nodep[m - 1], p);
}  /* buildnewtip */


void buildsimpletree(tree& curtree, const t_matrix_f& io, const vector<unsigned>& enterorder, unsigned nextsp)
{
	/* make and initialize a three-species tree */
	curtree.start=curtree.nodep[enterorder[0] - 1];
	setuptipf(curtree, io, enterorder[0], &curtree);
	setuptipf(curtree, io, enterorder[1], &curtree);
	hookup(curtree.nodep[enterorder[0] - 1], curtree.nodep[enterorder[1] - 1]);
	buildnewtip(curtree, io, enterorder[2], nextsp);
	insert_(curtree, curtree.nodep[enterorder[2] - 1]->back, curtree.nodep[enterorder[0] - 1],
			false);
}  /* buildsimpletree */


void addtraverse(tree &curtree, node *p, node *q, bool contin, unsigned& numtrees,
				 bool *succeeded, node *addwhere)
{
	/* traverse through a tree, finding best place to add p */
	insert_(curtree, p, q, true);
	numtrees++;
	if ( evaluate(curtree) > (curtree.best_tree->likelihood +
							   epsilonf * fabs(curtree.best_tree->likelihood)) )
	{
		copy_(&curtree, curtree.best_tree);
		addwhere = q;
		(*succeeded)=true;
	}
	copy_(curtree.prev_tree, &curtree);
	if ( !q->tip && contin )
	{
		addtraverse(curtree, p, q->next->back, contin,numtrees,succeeded, addwhere);
		addtraverse(curtree, p, q->next->next->back, contin,numtrees,succeeded, addwhere);
	}
}  /* addtraverse */


void re_move(tree& curtree, node **p, node **q)
{
	/* re_move p and record in q where it was */
	*q = (*p)->next->back;
	hookup(*q, (*p)->next->next->back);
	(*p)->next->back = NULL;
	(*p)->next->next->back = NULL;
	update(curtree, *q);
	update(curtree, (*q)->back);
}  /* re_move */


void globrearrange(tree& curtree, unsigned& numtrees, bool* succeeded)
{
	/* does global rearrangements */
	int j,k,num_sibs,num_sibs2;
	node *where,*sib_ptr,*sib_ptr2;
	double oldbestyet = curtree.likelihood;
	int success = false;

	tree globtree(curtree);
	tree oldtree(curtree);


	copy_(&curtree,&globtree);
	copy_(&curtree,&oldtree);
	for (unsigned i = curtree.spp ; i < curtree.nonodes ; i++ )
	{
		num_sibs = count_sibs(curtree.nodep[i]);
		sib_ptr  = curtree.nodep[i];
		for ( j = 0 ; j <= num_sibs ; j++ )
		{
			re_move(curtree, &sib_ptr,&where);
			copy_(&curtree,curtree.prev_tree);

			if ( where->tip )
			{
				copy_(&oldtree,&curtree);
				copy_(&oldtree,curtree.best_tree);
				sib_ptr=sib_ptr->next;
				continue;
			}
			else num_sibs2 = count_sibs(where);
			sib_ptr2 = where;
			for ( k = 0 ; k < num_sibs2 ; k++ )
			{
				node *addwhere = 0;
				addtraverse(curtree, sib_ptr,sib_ptr2->back,true,numtrees,succeeded, addwhere);
				if ( addwhere && where != addwhere && where->back != addwhere
					 && curtree.best_tree->likelihood > globtree.likelihood )
				{
					copy_(curtree.best_tree,&globtree);
					success = true;
				}
				sib_ptr2 = sib_ptr2->next;
			}
			copy_(&oldtree,&curtree);
			copy_(&oldtree,curtree.best_tree);
			sib_ptr = sib_ptr->next;
		}
	}
	copy_(&globtree,&curtree);
	copy_(&globtree,curtree.best_tree);
	if ( success && globtree.likelihood > oldbestyet )
	{
		*succeeded = true;
	}
	else
	{
		*succeeded = false;
	}
}


void rearrange(tree& curtree, node *p, unsigned& numtrees, unsigned *nextsp, bool *succeeded)
{
	// ignored
	node *addwhere = 0;
	node *q, *r;
	if ( !p->tip && !p->back->tip )
	{
		r = p->next->next;
		re_move(curtree, &r, &q);
		copy_(&curtree, curtree.prev_tree);
		addtraverse(curtree, r, q->next->back, false, numtrees,succeeded, addwhere);
		addtraverse(curtree, r, q->next->next->back, false, numtrees,succeeded, addwhere);
		copy_(curtree.best_tree, &curtree);
	}
	if ( !p->tip )
	{
		rearrange(curtree, p->next->back, numtrees,nextsp,succeeded);
		rearrange(curtree, p->next->next->back, numtrees,nextsp,succeeded);
	}
}  /* rearrange */



void summarize(tree& curtree, const t_matrix_f& io, double& sum_lengths, double& avg_percnt_stddev)
{
	/* print out branch lengths etc. */

	sum_lengths =  -curtree.likelihood;

	if (curtree.power == 2.0 && !curtree.minev )
	{
		unsigned missing_val = 0;
		for (unsigned i = 1; i < curtree.spp; ++i)
			for (unsigned j = 0; j < i; ++j)
				missing_val += io.is_missing(i, j);
		unsigned totalnum = curtree.spp * (curtree.spp - 1) - missing_val * 2;
		avg_percnt_stddev = 100.0 * sqrt(-curtree.likelihood / (totalnum - 2));
	}
}  /* summarize */



void nodeinit(node *p, bool lengths)
{
	/* initialize a node */

	for (unsigned i = 1; i <= 3; i++ )
	{
		std::fill(p->w.begin(),  p->w.end(), 1.0);
		std::fill(p->d.begin(),  p->d.end(), 0.0);
		p = p->next;
	}
	if ( (!lengths) || p->iter )
		p->v = 1.0;
	if ( (!lengths) || p->back->iter )
		p->back->v = 1.0;
}  /* nodeinit */


void initrav(node *p, bool lengths)
{
	/* traverse to initialize */
	if ( p->tip )
		return;
	nodeinit(p, lengths);
	initrav(p->next->back, lengths);
	initrav(p->next->next->back, lengths);
}  /* initrav */

void treevaluate(tree& curtree, const t_matrix_f& io, bool lengths)
{
	/* evaluate user-defined tree, iterating branch lengths */
	double oldlike;

	for (unsigned i = 1; i <= curtree.spp; i++ )
		setuptipf(curtree, io, i, &curtree);
	//	unroot(&curtree,nonodes2);							//OLDNEW

	initrav(curtree.start, lengths);
	if ( curtree.start->back != NULL )
	{
		initrav(curtree.start->back, lengths);
		evaluate(curtree);
		do
		{
			oldlike = curtree.likelihood;
			smooth(curtree, curtree.start);
			evaluate(curtree);
		} while ( fabs(curtree.likelihood - oldlike) > delta );
	}
	evaluate(curtree);
}  /* treevaluate */


void randomize_order(vector<unsigned>& enterorder)
{ /* randomize input order of species */
	std::reverse(enterorder.begin(), enterorder.end());
	unsigned spp = enterorder.size();
	for ( unsigned j = 0; j < spp; ++j )
		std::swap(enterorder[j], enterorder[rand() % spp]);
}



void maketree_usertree(tree& curtree, const t_matrix_f& io, istream& treefile, unsigned outgrno,
							const std::vector<std::string>& nayme, bool lengths)
{
	/* contruct the tree */

	// setup distances and weights
	for (unsigned i = 1; i <= curtree.spp; ++i )
		setuptipf(curtree, io, i, &curtree);

	treeread2 (treefile, &curtree.start, &curtree.nodep[0],
				   lengths, curtree.spp, nayme);
	curtree.start = curtree.nodep[outgrno]->back;
	treevaluate(curtree, io, lengths);
	//			clear_connections(&curtree,nonodes2);	//OLDNEW
}

bool maketree(tree& curtree, const t_matrix_f& io, unsigned& numtrees, unsigned spp,
				bool first_iteration, bool global, bool jumble,
				tree& besttree_over_all_runs)
{
	bool succeeded=false;

	//		nonodes2--;									//OLDNEW

	vector<unsigned> enterorder;
	for (unsigned i = 0; i < spp; i++ )
		enterorder.push_back(i + 1);
	if ( jumble )
		randomize_order(enterorder);


	unsigned nextsp = 3;
	buildsimpletree(curtree, io, enterorder, nextsp);
	curtree.start = curtree.nodep[enterorder[0] - 1]->back;
	nextsp = 4;
	while ( nextsp <= spp )
	{
		buildnewtip(curtree, io, enterorder[nextsp - 1], nextsp);
		copy_(&curtree, curtree.prev_tree);
		curtree.best_tree->likelihood = -99999.0;
		curtree.start = curtree.nodep[enterorder[0] - 1]->back;
		node *addwhere = 0;
		addtraverse(curtree, curtree.nodep[enterorder[nextsp - 1] - 1]->back,
					curtree.start, true, numtrees, &succeeded, addwhere);
		copy_(curtree.best_tree, &curtree);
		succeeded = true;
		while ( succeeded )
		{
			succeeded = false;
			curtree.start = curtree.nodep[enterorder[0] - 1]->back;
			if ( nextsp == spp  && global )
				globrearrange (curtree, numtrees, &succeeded);
			else
				rearrange(curtree, curtree.start,numtrees,&nextsp,&succeeded);
		}
		nextsp++;
	}

	// save as besttree_over_all_runs
	if (first_iteration ||
		besttree_over_all_runs.likelihood < curtree.best_tree->likelihood )
	{
		copy_(curtree.best_tree, &besttree_over_all_runs);
		return true;
	}
	return false;
}  /* maketree */



void fitch_make_nhstr(const t_matrix_f& io, t_fitch_params& params, t_fitch_results& results,
						t_progress_indicator& dots)
{
	unsigned spp = io.size();
	unsigned nonodes2 = spp * 2 - 2;

	assert(params.outgrno < io.size() );
	if (params.outgrno >= io.size() )
		params.outgrno = 0;

	// allocate space for nodes as well as leaves
	t_matrix_f node_matrix;

	node_matrix.matrix.allocate(nonodes2, 0.0);
	node_matrix.matrix.copy_from(io.matrix);


	tree curtree(spp, params.minev, params.negallowed, params.power, dots);

	tree besttree_over_all_runs(curtree);
	tree priortree(curtree);
	tree bestree(curtree);
	curtree.prev_tree = &priortree;
	curtree.best_tree = &bestree;



	unsigned numtrees = 1;

	for ( unsigned jumb = 0; jumb < params.njumble; jumb++ )
	{
		if (maketree(curtree, node_matrix, numtrees, spp, !jumb,
					params.global, params.jumble, besttree_over_all_runs))
			results.best_iteration = jumb;
		std::string indicator_str;
		if (curtree.dots.any_dots_printed())
			curtree.dots.print_indicator(ultoa(jumb + 1, indicator_str));
	}
	results.cnt_trees_searched = numtrees;

	// summarise
	copy_(&besttree_over_all_runs, &curtree);
	curtree.start = curtree.nodep[params.outgrno]->back;

	// printree(cout, curtree.start, false, curtree.spp, node_matrix.names);
	results.sum_lengths = 0.0;
	results.avg_percent_stddev = 0.0;
	summarize(curtree, node_matrix, results.sum_lengths, results.avg_percent_stddev);

	ostringstream ostrm;
	tree_to_nh_strm(ostrm, curtree.start, curtree.start, io.names);
	results.nh_tree = ostrm.str();
}


void optimize_nhstr(const string& tree_in_nh, const t_matrix_f& io,
					t_fitch_params_optimise& params, t_fitch_results& results,
					t_progress_indicator& dots)
{
	unsigned spp = io.size();
	unsigned nonodes2 = spp * 2 - 2;

	assert(params.outgrno < io.size() );
	if (params.outgrno >= io.size() )
		params.outgrno = 0;

	// allocate space for nodes as well as leaves
	t_matrix_f node_matrix;

	node_matrix.matrix.allocate(nonodes2, 0.0);
	node_matrix.matrix.copy_from(io.matrix);


	tree curtree(spp, params.minev, params.negallowed, params.power, dots);


	istringstream ifs(tree_in_nh);
   	maketree_usertree(curtree, node_matrix, ifs, params.outgrno, io.names, params.lengths);

	// printree(cout, curtree.start, false, curtree.spp, node_matrix.names);
	results.sum_lengths = 0.0;
	results.avg_percent_stddev = 0.0;
	summarize(curtree, node_matrix, results.sum_lengths, results.avg_percent_stddev);

	ostringstream ostrm;
	tree_to_nh_strm(ostrm, curtree.start, curtree.start, io.names);
	results.nh_tree = ostrm.str();
}



}
